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The Verwey phase transition in magnetite has been analyzed using the group theory methods. It 
is found that two order parameters with the symmetries X3 and A5 induce the structural transfor- 
mation from the high-temperature cubic to the low-temperature monoclinic phase. The coupling 
between the order parameters is described by the Landau free energy functional. The electronic 
and crystal structure for the cubic and monoclinic phases were optimized using the ab initio density 
functional method. The electronic structure calculations were performed within the generalized 
gradient approximation including the on-site interactions between 3d electrons at iron ions — the 
Coulomb element U and Hund's exchange J. Only when these local interactions are taken into ac- 
count, the phonon dispersion curves, obtained by the direct method for the cubic phase, reproduce 
the experimental data. It is shown that the interplay of local electron interations and the coupling 
to the lattice drives the phonon order parameters and is responsible for the opening of the gap at 
the Fermi energy. Thus, it is found that the metal-insulator transition in magnetite is promoted by 
local electron interactions, which significantly amplify the electron-phonon interaction and stabilize 
weak charge order coexisting with orbital order of the occupied tig states at Fe ions. This provides a 
scenario to understand the fundamental problem of the origin of the Verwey transition in magnetite. 

PACS numbers: 71.30.+h, 71.38.-k, 64.70.Kb, 75.50.Gg 



I. INTRODUCTION 

For many decades, magnetite — the oldest known mag- 
netic material — has been a fascinating subject of ex- 
tensive research. Besides its distinct magnetic proper- 
ties, the main interest was focused on the mechanism 
and physical consequences of the first-order phase tran- 
sition at Ty = 122 K, called the Verwey transition 
(VT), in which conductivity changes by about two orders 
of magnituteji For a long time it has been considered 
as an example of charge localization driven by metal- 
insulator transition^ in which ionic interactions deter- 
mine the electronic properties.- Only in the last decade 
the progress in experimental and theoretical techniques, 
as well as improved quality of samples, allowed one for a 
deeper insight into the nature of this transition^ Recent 
studies demonstrated that the VT is a cooperative phe- 
nomenon, in which an interplay between lattice, charge, 
and orbital degrees of freedom plays a decisive role. In 
this context, magnetite remains to be one of the most 
interesting materials among transition metal oxides^ 

At room temperature, magnetite (Fe304) crystalizes in 
the inverse spinel cubic structure with Fe ions occupying 
the tetrahedral A sites and octahedral B sites, as shown 
in Fig. [T] The ionic structure of magnetite in the inverse 
spinel phase is Fe 3+ (A)Fe 2 - 5+ (B)Fe 2 - 5+ (B)04". Below 
the Neel temperature Tjy = 851 K, magnetic moments at 
Fe ions are aligned antiparallely between A and B sites, 
which results in the ferrimagnetic state. Verwey pro- 
posed that the metal-insulator transition is caused by the 
charge order (CO) of Fe 2+ and Fe 3+ ions at the B sites, 
in planes perpendicular to the c axis.— The structural 
analysis, however, revealed the crystal distortion, which 
is incompatible with the Verwey model. In particular, 



the observation of half-integer reflections at (h, k, I + |) 
points indicated the doubling of the unit cell along the 
c axis.- Further diffraction studies by neutron* 7 ^ x-ray^ 
and electron^ scattering established the monoclinic sym- 
metry of the low temperature (LT) phase realized below 
Ty ■ Indeed, the number of inequivalent A and B atomic 
positions found by nuclear magnetic resonance (NMR) 
for the phase below the V T 11 ' 12 ' 13 agrees with the mon- 
oclinic structure. 

The models to describe the CO in magnetite, consis- 
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FIG. 1: (Color online) The crystal structure of magnetite in 
the cubic FdZm symmetry. Iron and oxygen ions are repre- 
sented by grey and dark (orange and blue) balls. Fe(_B) ions 
have six O neighbors, while each O ion has three Fe(B) and 
one Fe(A) neighbor. 
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tent with the LT phase of Fe304, have been proposed in 
many studies i 10 ' 11 ' 12 ' 13 High-resolution neutron and x- 
ray scattering measurement o 14 ' 15 revealed the distribu- 
tion of Fe-0 bond distances in the BO§ octahedra (be- 
tween 1.96 and 2.11 A), being the indication of charge 
disproportionation on Fe(B) ions. The CO proposed by 
Attfield et al^>^- is considerably more complex than 
the Verwey modeli itself, and consists of four different 
B sites in a unit cell, with approximately two valence 
states: Fe +24 and Fe +2 6 . Providing a convincing evi- 
dence in favor of such a state has been a challenge for 
several direct probing methods including the resonant x- 
ray scattering. Despite preliminary negative results ; 16 ' 17 
the evidence that the CO is indeed fractional has accu- 
mulated recently . 18 ' 19 ' 20 ' 21 It is clear from the character 
of the observed CO that a simple ionic mechanism sug- 
gested early on by Anderson 3 - cannot explain the VT. 

The electronic structure of magnetite was studied 
by the density functional theory (DFT) calculations in 
the local spin density approximation (LSD A).— This 
approach is based on the local density approximation 
(LDA) and provides a basic description of the electronic 
structure and magnetic properties of Fe304 in the cu- 
bic phase above Ty . 22 ' 23 ' 24 The CO instability induced 
by local on-site and intersite Coulomb interactions was 
studied in the cubic symmetry within the LSDA+J7 
method . 25 ' 26 For a distorted (orthorhombic) symmetry, 
the self-interaction correctional approach was used, and 
it has been shown that the CO proposed by Verwey does 
not appear in the ground state of magnetite. The calcu- 
lations for the experimentally observed monoclinic struc- 
ture, performed within the LSDA+J7 (see Refs. HHHl) 
and generalized gradient approximation (GGA) with lo- 
cal Coulomb interaction U (GGA+[7, see, e.g. Ref. [30h 
approaches proved the existence of fractional charge dis- 
proportionation in fair agreement with the experimental 
analysis fi 4 In addition, these studies revealed the pres- 
ence of the orbital order (00) of the occupied t 2g states 
at .B-sites t 29 ' 30 in agreement with recent experiments per- 
formed on thin films of Fe304J 3 -^ The obtained magnitude 
of the insulating gap agrees well with the photoemision 
data below 7V , 32 ' 33 ' 34 On the contrary, in absence of elec- 
tron interactions (for U — 0) the electronic state remains 
metalic even in the distorted structure, showing the es- 
sential role played by local electron interactions in the 
VT. These ab initio calculations demonstrate also that 
the microscopic understanding of the VT requires to in- 
volve lattice distortions and electron-phonon (EP) inter- 
action. 

The experimental evidence that phonons participate in 
the mechanism of the VT is very convincing. The struc- 
tural analysis showed that the LT phase results from the 
condensation of modes at three reciprocal lattice points: 
k r = (0,0,0), k A = (0,0, i), and k x = (0,0, l)^ 
in units of — , where a is lattice constant of the cubic 
structure. The pronounced softening of the C44 elas- 
tic constant ^ observed from room temperature down 
to Ty, has been explained by the coupling between the 



shear mode and the charge-density T 2g order parame- 
ter (OP) at the r point i^ 6 - The decrease of C44 is asso- 
ciated with the softening of the surface Rayleigh mode 
observed by the Brillouin scattering^ In the bulk, no 
phonon softening has been observed^ however, critical 
fluctuations revealed by neutron scattering strongly in- 
dicate the existence of phonon precursor effects which 
drive the structural deformation. First, the spotlike dif- 
fuse signal was observed at the reciprocal space point 
kA = (4, 0, i)^ 9 - Yamada proposed that the VT re- 
sults from the condensation of a coupled charge density- 
phonon mode with A5 symmetry, explaining the intensi- 
ties of neutron signal and kA Bragg reflections Second, 
a new type of planar diffuse scattering was found over a 
broad region of the reciprocal space at high temperature 
T v <T <T V + 100 K^iii 2 . The signal around the T and 
equivalent points was interpreted as the Huang scattering 
due to the local strain field, and was studied in the molec- 
ular polaron modeL 43 ' 44 The analysis of intensities at kx 
points led to the conclusion that transverse phonons with 
the X 3 symmetry dominate in neutron scattering J 5 - The 
phonon mechanism is further supported by oxygen iso- 
tope effect^ 6 - Raman ; 47 ' 48 ' 49 the extended x-ray absorp- 
tion fine structure (EXAFS)^ and by nuclear inelastic 
scattering (NIS) measurements.— 

The idea that both Coulomb interactions and the EP 
coupling are responsible for the VT was first explored 
by Ihle and LorentziS 2 - They combined the tight-binding 
model introduced by Cullen and Callen 53 with the Ya- 
mada model^ In their approach, the CO below Ty is 
driven by the intersite electrostatic interactions and the 
condensation of the A5 phonon^ 4 - Above Ty, the model 
explains the short-range polaronic order and provides a 
good description of the electric conductivity^ 5 - Indeed, 
the polaronic character of charge carriers has been ob- 
served in the optical conductivit y 48 ' 56 and in photoemi- 
sion studies* 3 ^ The cooperative mechanism of the VT was 
next studied also in the Peierls-Hubbard model; 57 where 
strong on-site interaction induces the OO and lattice dis- 
tortion. Since the active electrons on Fe +2 ions occupy 
partly filled degenerate t 2g states, the EP coupling could 
be further enhanced by the Jahn- Teller effect. This mech- 
anism was confirmed by the dynamical mean-field theory 
(DMFT) combined with the DFT in the GGA+DMFT 
method^ and by GGA+C/ calculations, 59 which showed 
that the local Jahn- Teller distortions modify the elec- 
tronic structure lead indeed to opening of the insulating 
gap. 

In the previous paper p! - we reported the phonon spec- 
trum obtained using the direct method? 6 - 1 - and analyzed 
the EP interactions in presence of local Coulomb interac- 
tion U. We have found that two phonon modes with X3 
and A5 symmetry, which significantly distort the 5-sitc 
octaedra, couple strongly to t 2g electronic states. The 
X3 phonon opens the gap at the Fermi energy, and drives 
the metal-insulator transition. Using the group theory, 
we proved that these two modes are primary order pa- 
rameters describing the structural phase transition. 
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The purpose of the present paper is to focus on the 
details of the group theory analysis of the phase transfor- 
mation in the magnetite, from the cubic FdSm to mon- 
oclinic P2/c symmetry, and to relate it to the simulta- 
neous changes in the electronic structure and lattice dy- 
namics. A full list of relevant OPs will be presented and 
their couplings will be discussed on the basis of Landau 
free energy. In the second part of the paper we focus on 
the results of structural optimization and lattice dynam- 
ics calculations in presence of local electron interactions, 
and demonstrate that the phase transformation in the 
VT involves the EP interaction. 

The paper is organized as follows. In Sec. II the group 
theory analysis of the VT is presented. Sec. Ill is de- 
voted to the electronic structure calculation, where we 
first present the method used (Sec. IIII A|) . and next an- 
alyze the optimized crystal structures (Sec. IIIIB|) . In 
Sec. IV we describe lattice dynamics and classify the 
obtained modes for the magnetite using the irreducible 
representations of the cubic group. These results serve 
as a basis to discuss the EP interactions and their effect 
on the electronic structure in Sec. V. We show a strong 
coupling between the electronic properties and lattice dy- 
namics, and in this context discuss the mechanism of the 
VT in Sec. VI. The paper is concluded in Sec. VII, where 
we provide a coherent view on the VT in magnetite and 
summarize the presented results. 



II. GROUP THEORY ANALYSIS 

A. Primary and secondary order parameters 

In a structural phase transition the symmetry is low- 
ered from the high-symmetry phase, usually realized at 
high temperature, to the low-symmetry one, stable at 
low temperature. The Verwey transition in magnetite 
belongs to this class of phase transitions, so we consider 
first the Landau free energy functional to characterize the 
OPs responsible for the observed lowering of symmetry. 
In the Landau-type phase transition the space group £ of 
the low-symmetry phase is a subgroup of the space group 
H of the high-symmetry phase, i.e., £ C H. There ex- 
ists often one irreducible representation (IR) of 7i, which 
reduces the symmetry to £. This allows one to identify 
an OP as a vector in an IR space, T^jfa), where k is a 
wave vector of the irreducible star, j is the index of ray 
representation, and r\i are the components of the OP. 62 
The OP which defines the transition and determines the 
symmetry of the low-symmetry phase is called a 'primary 
OP, and the symmetry reduction diagram for it can be 
written as follows 

w-tTkjfai)]-/:. (i) 

In contrast, a secondary OP yk,j{T)i) is associated with 
an IR which reduces the symmetry of the high-symmetry 
space group to an intermediate space group I, such that 



X becomes a subgroup to Tl, and a supergroup for £, i.e., 
£ C T C TLf& Thus, the symmetry reduction diagram 
has the form 

[jkAvi)} ->i. (2) 

The secondary OP appears in the Landau free energy ex- 
pansion always in a term which provides a linear coupling 
of this secondary OP to a quadratic (or higher power) ex- 
pression in the relevant primary OP. At the phase tran- 
sition the primary OP starts to develop, while the sec- 
ondary OPs also become finite because of their coupling 
to the primary OP. 

In principle a phase transition could be induced by 
two or more primary OPs, and they can be coupled 
to the secondary ones. Thus, in general N primary 
OPs, IRi,IR2,- • • ,IR n, reduce the high-symmetry space 
group H into £i, £2, £n low-symmetry space 

groups, respectively, 

H -> [r kl ,j, (m)] -> A, 
ft [ r k 2 j 2 (^)] -> £2, 

n -[Tk^fo)]-/:*. (3) 

Due to a non- linear coupling between the primary OPs, 
the final low-symmetry space group £ is an intersection 
of all low-symmetry space groups 

£ = £ 1 n£ 2 r\...£ N . (4) 

Thus, the low-symmetry space group £ consists of those 
symmetry elements which are present in all subgroups 
£1, £2, ■ ■ ■ ,£n- The secondary OPs are associated with 
IRs which reduce the symmetry of the high-symmetry 
space group to intermediate space groups T m , such that 



TABLE I: List of OPs from the parent space group Fd3m 
(No=227) and basis {(1, 0, 0), (0, 1, 0), (0, 0, 1)} to the mon- 
oclinic phase P2/c (No=13, unique axis b, choice 1), basis 
{(±,-±,0),(±,±,0), (0,0,2)} and origin (±,0,±) relative to 
the original face center cubic lattice. Size is the ratio of the 
volumes of primitive low-symmetry to high-symmetry unit 
cells. The direction of the IR vector given in the last column 
indicates its nonvanishing components as well as the relations 
between the components of the OPs. 



IR 


size 


subgroup 


No 


IR vector 


rt 


1 


Fd3m 


227 


(a) 


it 


1 


IA\/amd 


141 


(6,0) 


it 


1 


C2/m 


12 


(c, c,0) 


it 


1 


Imma 


74 


(d,0,0) 


it 


1 


C2/m 


12 


(d, e, — e) 


Xi 


2 


Pmma 


51 


(/, 0,0, 0,0,0) 


x 3 


2 


Pmna 


53 


(3,0,0,0,0,0) 
(0, 0,0,0, h, -ifc) 


A 2 


4 


Pcca 


54 


A 4 


4 


Pcca 


54 


(0,0,0,0,p,MP) 


A 5 


4 


Pbcm 


57 


(0, 0, 0, 0, 0, 0, 0, 0, q, —nq, -fiq, -q) 
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TABLE II: Subgroup relationships and the IRs creating the symmetry reduction from space group Ti to a given T and £ space 
groups, respectively. In brackets are the sizes of primitive unit cells. 



Subgroup relationships Group No's Secondary OP Primary OP 

Pd3m(l) 5 Mjomdfl) D Pmma(2) 5 P6cm(4) 227, 141, 51, 57 r+, Xi Aj 

Fd3m{l) D /mma(l) D Pmma(2) D Pcm(4) 227, 74, 51, 54 T+, Xr A 2 , A 4 

Fd3m(l) D J4iamd(l) D Pmna(2) 227, 141, 53 P^ X 3 

Fd3m(l) D Jmmo(l) D Pmna(2) 227,74,53 r+(d,0,0) X 3 

Fd3m(l) D Jmma(l) D C2/m(l) 227,74,12 r+(d,0,0) T+ (d,e,-e) 



TABLE III: 


Intersections of subg 


roups listed in Table I which lead to low-symmetry space 


group P2/c. In all cases the size is 


equal 4, the basis vectors are {(| 


-i,0),(i,i,0), (0,0,2)} and the origin is (|,0,±). 




No 


IRs 


subgroups giving P2/c 


subgroup No's 


1 


X 3 ,A 5 


Pmna n Pbcm 


(53, 57) = 13 


2 


X 3 ,A 2 


Pmna n Pcca 


(53, 54) = 13 


3 


X 3 ,A 4 


Pmna n Pcca 


(53, 54) = 13 


4 


x 3 ,A 5 ,r+ 


Pmna n Pbcm n C2/m 


(53, 57, 12)= 13 


5 


x 3 ,A 2 ,r+ 


Pmna n Pcco n C2/m 


(53, 54, 12)= 13 


6 


x 3 ,A 4 ,r+ 


Pmna n Pcca n C2/m 


(53, 54, 12)= 13 


7 


X 3 ,A 5 ,A 2 


Pmna n Pbcm n Pcca 


(53, 57, 54)= 13 


8 


X 3 ,A 5 ,A 4 


Pmna n Pbcm n Pcca 


(53, 57, 54)= 13 


9 


X 3 ,A 2 ,A 4 


Pmna n Pcca n Pcca 


(53, 54, 54)= 13 


10 


X 3 ,A 5 ,A 2 ,P 5 + 


Pmna n Pbcm n Pcca n C2/m 


(53, 57, 54, 12) = 13 


11 


X 3 ,A 5 ,A 4 ,r+ 


Pmna n Pbcm n Pcca n C2/m 


(53, 57, 54, 12) = 13 


12 


X 3 ,A 2 ,A 4 ,r+ 


Pmna n Pcca n Pcca n C2/m 


(53, 54, 54, 12) = 13 


13 


X 3 ,A 5 ,A 2 ,A 4 


Pmna n Pbcm n Pcca n Pcca 


(53, 57, 54, 54) = 13 


14 


X 3 ,A 5 ,A 2 ,A 4 ,r+ 


Pmna n Pbcm n Pcca n Pcca nC2/m 


(53, 57, 54, 54, 12)= 13 



L n C 1 m C H, where n — 1, 2, . . . N and m should be 
found by inspection the group-subgroup relationships 

ft -> bk u n(Vi)] 

ft -> hfc 2 ,h(Vi)] ~>Z2, 



H 



[7k„ 



Xm)] 



(•5) 



For instance, a symmetry reduction in shape memory al- 
loy NiTi, associated with the symmetry change PtoSto —> 
P2\/m 1 is a phase transition which involves two primary 



OP; 
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As we demonstrate in this Section, the sym- 



metry reduction realized in the VT of magnetite requires 
also multiple primary and secondary OPs. 

The high-symmetry space group of Fe304 above Ty is 
H = Fd3m (Fig. 1). The low-symmetry phase £ was 
identified as the monoclinic Cc symmetryi£&i£ with a 
unit cell of approximate size y/2a x \/2a x 2a. As argued 
by Iizumi et ai.^- a smaller cell aj\[2 x a/V2 x 2a with 
the Pmca symmetry can be used for refinement almost 
all Bragg reflections. In Ref. [HI, the monoclinic sym- 
metry P2/c with orthorhombic constraints was used in 
the refinement procedure. This symmetry describes all 
Bragg peaks except for three very weak reflections arising 
from atomic displacements smaller than 0.01 Aji£ Indeed, 
the unit cell used in the refinement procedure takes into 
account all atomic displacements coming from kr, kA, 
and kx points. Therefore, the P2/c symmetry explains 



majority of physical properties connected with the VT. 

Here, we present a detailed group theory analysis of 
the structural transition for C = P2/c (No=13), with ba- 
sis {(I, -I, 0), (A, ±,0), (0,0, 2)} and the origin (±,0,±) 
relative to the original face center cubic lattice. Us- 
ing the COPL and ISOTROPY computer code a 66 i 67 we ob- 
tained the lists of possible intermediate {T m } and low- 
symmetry {£«} space groups, and show them in Table I. 
The space groups, which have been found, are reduced by 
IRs from the high-symmetry space group TL = FdZm. As 
is clearly seen, there is no IR which reduces Fd3m to low- 
symmetry C — P2/c. This means that the phase transi- 
tion in Fe3C>4 is driven by at least two primary OPs. A 
closer inspection of space group-subgroup relationships, 
partly shown in Table II, leads us to the conclusion that 
in principle there are five primary OPs: X3, A2, A4, and 
A5, and r^" (T2 9 ) within direction (d, e, — e). Indeed, the 
intersection of space groups being the result of symmetry 
reduction of Fd3m by the above mentioned IRs provides 
the observed low-symmetry space group P2/c, Table III. 

In the the simplest case, the P2/c symmetry can be 
induced by two OPs. Taking the X3 and A5 symmetries, 
the reduction diagrams can be written 



Fd3m -> [X 3 , k = (0, 0, 1)] ~> Pmna, 



Fd3m 



A,,k 



0,0, 



Pbcm. 



(6) 



Now, if we take common symmetry elements (an inter- 
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section) of these two groups, we get 

Pmna n Pbcm = P2/c, (7) 

so the low-symmetry space group has indeed P2/c sym- 
metry. The remaining IRs involved in the phase transi- 

tion r+ (A lg ), r+ {E g ), r+ (r lg ), r+ (r 2g ) within di- 

rection (d, 0, 0), and X\ are the secondary OPs. On the 
one hand, it is evident that can be either primary or 
secondary OP. If only one component d is nonzero it has 
a secondary character and its reduction diagram reads 

Fd3m^ [r+(<i,0,0),k= (0,0,0)] -> Imma. (8) 

The resulting group fulfills the following subgroup rela- 
tionship 

Fd3m D Imma D Pmna, (9) 

so Tg" as the secondary OP couples to the primary OP 
A3. On the other hand, if all three components of 
are nonzero, it has primary character with the reduction 
diagram 

Fdlm -> [T+(d, e, -e), k = (0, 0, 0)] -> C2/m. (10) 

Because of the following relationship 

Fdim D Imma D C2/m, (11) 

r^"(<i, 0, 0) is the secondary parameter of the T§(d, e, — e) 
OP. In what follows, we shall use the point group notation 
of the IRs at the T point. This notation is usually used 
in experimental works on phonon spectroscopy. 

As is clear from this analysis, the A3 OP plays an 
exceptional role among other symmetries, being present 
in all intersections listed in Table III. So one recognizes 
that this OP is essential to generate the LT monoclinic 
phase and it confirms its precursor behavior found in the 
diffuse scattering study. As we showed by the ab initio 
calculations, this mode strongly couples to electrons and 
induces the metal-insulator transition^ 

At the kA point, there are three possible OPs: A2, A4, 
and A5. In principle, by coupling to A3 either of these 
modes can induce P2/c symmetry, but one of them, A5, 
has a strong support from the experiment. Originally 
proposed by Yamada^S — this mode was found to be the 
main component of lattice distortion in critical neutron 
scattering *2£ Neutron studies by Iizumi et al& confirmed 
the A5 pattern of displacements, which gives rise to kA 
Bragg peaks. The EP interaction for the A5 mode was 
studied within the model of Ihle and Lorenz^ At the 
r point there could be only one primary OP with T2 g 
symmetry. It is related with the softening of the elastic 
constant C44. Notice that in the cubic system one has 
C44 = C55 = C66, and therefore one is unable to differenti- 
ate between the secondary [d, 0, 0) or primary (d,e,—e) 



OP characters. A strong EP interaction for the Ti g mode 
was found by the Raman experiments 9 , and in neutron 
diffuse measurements.— Taking into account that only 
t2 g electronic states are present near the Fermi level, a 
strong coupling to phonons with the same symmetry is 
not surprising. Altogether, from experimental observa- 
tions and from ab initio calculations, presented in detail 
in the next Sections, we have strong evidence for the fol- 
lowing three primary OPs being involved in the phase 
transition: A 3 , A 5 , and T 2g . 

To appreciate fully the complexity of the resulting LT 
phase, we emphasize that although the LT symmetry 
is uniquely determined by the primary OPs, the exact 
atomic positions result from all lattice distortions, asso- 
ciated also with secondary OPs. The X\ type distortion 
was reported in x-ray^ and in neutron^ diffraction. The 
mechanism stabilizing the incommensurate mode A5 at 
wavevector kA by simultanuous condensation with the 
X\ phonon was discussed by Iizumii^ 9 - The X\ mode in- 
duces the longitudinal displacements with X — a observed 
in the monoclinic structured There are a few secondary 
OPs at the zone center, and they are likely to partic- 
ipate in the LT structure as welled For example, the 
anomaly related to the observed changes in frequency 
and linewidth at the phase transition was found for the 
highest Raman A\ g modeM- 



B. Free energy functional 

Knowing the primary and secondary OPs and their 
IRs, one can construct the invariant quantities which 
do not change under symmetry group transformations, 
and next expand the free energy T into a series of the 
components which involve different OPs^i Each compo- 
nent is a linear combination of its electronic and phononic 
parts, which represent the electronic density and atomic 
displacement contributions, respectively. For instance, 
the component of A3 can be written in the form which 
includes its electronic and phononic contribution: g — 
9c\ + .9ph- Note that in the present case of active ti g de- 
grees of freedom of 3d electrons at iron ions, the electronic 
part g e \ may as well involve the orbital occupations. 

We limit the expansion of the free energy T to the 4-th 
order terms, with the exception of two primary OPs A3 
and A 5 , for which we include also terms of 6-th order as 
they are typical for the discontinuous (first-order) phase 
transitions. A complete expansion gives too lengthy ex- 
presion to be reproduced here, so we leave in the expan- 
sion only components of OPs required by symmetry as 
shown in Table I, and write down only the potentially 
nonvanishing terms. Each term of the series must be 
multiplied by factors depending on thermodynamical pa- 
rameters, however, we shall not discuss this aspect here. 
The free energy reads 
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T = (g 2 + q 2 + g 4 + q 4 + g 2 q 2 + g 6 + .gV + <?V + g 6 ) + (p 2 + h 2 + g 2 P 2 + g 2 h 2 + q 2 P 2 + q 2 h 2 + P 2 h 2 + gq (h + p)) 

+ if + fq 2 + fp 2 + fh 2 + fph + fg 2 ) + (a 2 + ag 2 + aq 2 + ap 2 + ah 2 ) 

+ (b 2 + b 3 + bg 2 + bq 2 + bp 2 + bh 2 ) + (c 2 + c 2 g 2 + c 2 q 2 ) + (d 2 + e 2 + dg 2 + dp 2 + dh 2 + dph + ep 2 + eh 2 ). (12) 



First bracket in Eq. (fT2|) corresponds to the X3 and 
A5 IRs and their couplings (see Table I). We anticipate 
that the lowest coupling term between them g 2 q 2 must be 
large and should lower the free energy T . Second bracket 
correponds to A2 and A4 IRs and to the coupling terms 
among them. The second line in Eq. (fl2|) describes the 
coupling of IR X\ with the primary OPs, being of linear- 
quadratic form. The interaction between the X\ and A5 
phonons discussed by Iizumi^ 9 - is described by the fq 2 
term. The following terms given in fourth, fiveth, sixth 
and seventh brackets determine the coupling of IRs of the 
crystal point groups A± g , E g , T\ g , and T2 g , respectively, 
with the primary OPs. Again their forms are adequate to 
describe the influence of secondary OPs on the transition. 
In particular, the term dg 2 in the last line, describing the 
coupling between the X 3 and T 2g OPs, contributes to the 
softening of the C44 elastic constant. 3 - Other couplings 
involving the zone-center IRs may account for anomalies 
observed by the Raman studies * 48 ! 49 

In order to identify the relevant part of the free energy 
Eq. (|12p which drives the VT in magnetite, we consider 
next the changes of the electronic structure (Sec. Ill) 
and the phonon spectra (Sec. IV) at the transition. In 
fact, we shall demonstrate below that these properties 
are strongly coupled to each other, and the VT is driven 
by the electron-phonon interaction enhanced by the local 
Coulomb interaction U, see Sec. V. 



III. ELECTRONIC STRUCTURE 

A. Computational methods 

The lattice parameters of Fe304 were optimized to- 
gether with the electronic structure using the total en- 
ergy DFT approach. The numerical calculations were 
performed by means of the VASP program^ within the 
GGA+U approach^ The program uses highly accurate 
full-potential projector-augmented wave (PAW) method, 
originally proposed by Blochl— and implemented by 
Kresse and Jouberti^ 3 - The wavefunctions in the core re- 
gion are obtained by the linear transformation from vari- 
ational pseudo- wavefunctions. In the augmented region 
the wavefunctions are expanded in the plane wave ba- 
sis. Two supercells with the periodic boundary condi- 
tions were chosen: ax ax a and a/v2 x a/y/2 x 2a for the 
FdZm and P2/c structures, respectively, with 56 atoms 
each. The latter supercell was used also for the reference 
system of FdZm symmetry, when total energies of two 
considered structures were compared. 



For a given ionic configuration, the Kohn-Sham Hamil- 
tonian was diagonalized by the iterative residual mini- 
mization method, - which allows one for efficient opti- 
mization of the charge density and wavefunctions. In 
the initial steps also the blocked Davidson minimization 
method was applied. To improve convergence a small 
smearing of electronic states at the Fermi surface was 
introduced with the effective parameter a — 0.2 eV. 
The Monkhorst-Pack scheme^ was applied for summa- 
tion over two k-point grids 6x6x6 and 6x6x2 relevant 
for Fd3m and P2/c symmetries, respectively, which are 
large enough for a meaningful energy comparison. In Rcf. 

we used the smaller grids, i.e., 4x4x4 and 4x4x2, 
which give already satisfactory values for lattice param- 
eters and energies. However, larger k-point grids were 
used at present to reach a higher accuracy as they give 
even more reliable values of force constants and phonon 
frequencies. The energy cut-off for the plane wave ex- 
pansion was set at 520 eV after we have verified that the 
contribution of states with higher energies was negligible. 

The local electron interations between 3d electrons on 
Fe ions were included in the Hartree-Fock approximation, 
as usually done in the LDA+J7 method^ As a result one 
finds the total energy of the system E to t in the form 

Etot — Eqqa + Ejj — Edc, (13) 

where £gga is the energy obtained in the GGA ap- 
proach, Ejj describes the contribution due to local in- 
teractions parametrized by the Coulomb element U and 
Hund's exchange J, and E^ c is the double counting cor- 
rection term, i.e., the averaged electron-electron inter- 
action, which has to be subtracted. The details of the 
Hamiltonian used in the calculations were presented in 
Ref. [zS 

In the present calculations we have chosen the follow- 
ing parameters: U = 4.0 eV and J = 0.8 eV. The value 
of U agrees with the constrained DFT calculations ) 23 ! 58 
and is somewhat reduced from the ionic value of 6.4 eV 
estimated by Zaanen and Sawatzky. 7 & The value of J 
was obtained from the atomic values of the Racah pa- 
rameters for Fe 2+ ions-^ B = 0.131 eV and C = 0.484 
eV. Here we use an average value of Hund's exchange,-^ 
J = -B + C, as usually implemented within the LDA+/7 
method^ 8 . Note that the value of J = 1.0 eV used in 
some LSDA+t/ calculations! 9 - and in the LDA+DMFT 
method 5 ^ is somewhat overestimated (in fact this value 
applies to a pair of e g electrons rather than to a pair of t 2g 
electrons). The nearest-neighbor Coulomb interaction is 
one order of magnitude smaller V ~ 0.3 — 0.4 eV^ and 
we assume that it is to a large extent included already in 
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the GGA scheme. 

The optimal atomic configuration, which gives the 
minimum of the total energy, E = min{E tot }, was 
found using the conjugate gradient and quasi-Newton 
procedures^ Because of phonon calculations, which are 
based on the well-optimized supercell, the terminating 
criteria for the electronic and ionic degrees of freedom 
were very strict: 10~ 7 eV and 10~ 5 eV, respectively. To 
get more precise atomic positions we used also force cri- 
terion, in quasi-Newton procedure, which allowed us to 
get residual forces less than 10 me V/A (the pressure 
was then less than 0.1 kbar). 

Phonon frequencies were calculated only for the cu- 
bic structure, using the direct method^ implemented 
in the PHONON program.— According to the Hellmann- 
Feynmann (HF) theorem, the atomic forces can be ob- 
tained by displacing atoms from their equilibrium posi- 
tions. The minimum number of displacements depends 
on crystal symmetry and on the number of nonequivalent 
atoms. For magnetite only three independent displace- 
ments for Fe(A), Fe(B), and O atoms are sufficient. To 
minimize systematic errors we performed displacements 
in positive and negative directions. 

The combination of derived displacements and forces 
allows one to obtain the force constants matrix elements 
by the singular value decomposition method. A spe- 
cial treatment was applied for calculations with finite U, 
where the convergence of the electronic part was in some 
cases very slow. When this occured, the calculations of 
the HF forces were initialized using the wavefunctions op- 
timized for the cubic structure. For magnetite, the direct 
method provides exact frequencies at the F and X points. 
The phonon frequencies at other points were evaluated 
with only small errors due to large supercell used in the 
present calculations. 



B. Crystal optimization 

The structural properties of the cubic and monoclinic 
phases have been discussed in many previous studies. In 
Tabs. IV and V we summarize the electronic and lattice 
parameters obtained in the present calculations, compar- 
ing them with the experimental data. We focused mainly 
on the effects connected with the Coulomb interaction U. 
For the cubic symmetry, the lattice constant a and the 
internal oxygen parameter x depend rather weakly on U, 
and show overall good agreement with experimental data 
of Refs. mmiHii (Table IV). As expected and noticed 
before^ the magnetic moments on particular atoms in- 
crease with U, improving the agreement between the elec- 
tronic structure calculations and the experimental data. 
The average magnetic moment per one formula unit m av , 
however, does not change with U and J. It shows that 
magnetic polarization in this material is already at maxi- 
mum and can be therefore well described in the band pic- 
ture, while reliable calculations concerning the magnetic 
polarization of individual Fe ions require local electron 



TABLE IV: Ground state lattice and electronic parameters 
for Fd3m and P2/c structures. The magnetic moments of 
individual ions Fe(A) (rriA), Fe(B) (tub), and O (mo) ob- 
tained within the present GGA (U = J = 0) and GGA+U ( 
with U = 4.0 eV and J = 0.8 eV) calculations for the high- 
temperature Fd3m phase are compared with the experimental 
data, wherever available. The lattice constants (a, b, c) are in 
A, magnetic moments fj,i in /is, and energies E to t and the gap 
A 9 in eV. 



phase 


quantity 


GGA 


GGA+f/ 


experiment 


Fd3m 


a 


8.377 


8.446 


8.394° 




X 


0.2545 


0.2549 


0.2549° 




m A 


-3.40 


-3.97 


-3.82 6 




TUB 


3.52 


3.87 






mo 


0.08 


0.04 






m av 


3.96 


3.96 


4.05 c 




-Etot 


-426.39 


-386.80 






A 9 










P2/c 


a 


5.926 


5.967 


5.944° 




b 


5.926 


5.995 


5.925° 




c 


16.752 


17.034 


16.775° 




P 


90.004 


90.417 


90.237° 




-Etot 


-426.63 


-389.28 






A 9 





0.33 


0.14 d 



° Reference 


15 




b Reference! 


81 




c Reference] 


82 




d Reference 


56 





interactions {U, J}. One finds that magnetite is metallic 
(A g = 0) in the cubic phase, independently of the ac- 
tual values of U and J. Note that the total ground state 
energy of the crystal -Etot increases with U, but direct 
comparison between the energies obtained in the GGA 
and the GGA+U is meaningless. 

The P2/c structure was optimized starting from the 
experimentally determined geometry,™ but we did not 
impose orthorhombic Pmca symmetry constraints on 
atomic positions. The theoretical crystal structures, ob- 
tained also by the VASP program, were already discussed 
in detail in Ref. HU. We note that some discrepancies 
between the results of these two calculations may follow 
from different energy cut-off and termination conditions, 
as well as from slightly different values of U and J used. 
In the GGA calculations performed with U = J = 0, the 
lattice constants, atomic positions, and the total energy 
converge to those of the high-symmetry cubic structure 
(Tabs. IV and V). This behavior agrees with the molecu- 
lar dynamics simulations presented in Ref. 59. The prox- 
imity of the total energies -Etot obtained in both phases 
indicates that by going beyond the GGA one of these 
structures will be stabilized. Note that the energy dif- 
ference between these two phases obtained in the GGA 
calculations is too small to arrive at definite conclusion 
concerning the stability of the LT phase of P2 /c symme- 
try. 

Indeed, the optimization performed with U = 4.0 eV 
and J — 0.8 eV gives lower total energy of the P2/c 
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TABLE V: Atomic positions for iron ions (at A and B sites) and oxygen ions in the P2/c structure, as obtained in the GGA 
calculations (left, U — J = 0) and in the GGA+U calculations (middle, with U = 4.0 eV and J = 0.8 eV), compared with the 
experimental data of Ref. [15 (right) . 







U = J = 




U = 


4.0 eV, J = Oi 


5 eV 




Experiment 




atom 


X 


y 


z 


X 


y 


z 


X 


y 


z 


Al 


0.2490 


0.0 


0.0625 


0.2491 


0.0104 


0.0637 


0.25 


0.0034 


0.0637 


A2 


0.25 


0.5 


0.1875 


0.2493 


0.4991 


0.1885 


0.25 


0.5061 


0.1887 


Bla 


0.0 


0.5 


0.0 


0.0 


0.5 


0.0 


0.0 


0.5 


0.0 


Bib 


0.5 


0.5 


0.0 


0.5 


0.5 


0.0 


0.5 


0.5 


0.0 


B2a 


0.0 


0.0 


0.25 


0.0 


0.0050 


0.25 


0.0 


0.0096 


0.25 


B2b 


0.5 


0.0 


0.25 


0.5 


0.9996 


0.25 


0.5 


0.0096 


0.25 


B3 


0.25 


0.25 


0.375 


0.2502 


0.2604 


0.3804 


0.25 


0.2659 


0.3801 


54 


0.25 


0.75 


0.375 


0.2529 


0.7553 


0.3740 


0.25 


0.7520 


0.3766 


01 


0.25 


0.2591 


-0.0023 


0.2508 


0.2689 


0.0000 


0.25 


0.2637 


-0.0023 


02 


0.25 


0.7408 


-0.0023 


0.2491 


0.7572 


-0.0026 


0.25 


0.7461 


-0.0029 


03 


0.25 


0.2409 


0.2523 


0.2487 


0.2370 


0.2541 


0.25 


0.2447 


0.2542 


04 


0.25 


0.7591 


0.2523 


0.2482 


0.7708 


0.2474 


0.25 


0.7738 


0.2525 


05a 


-0.0091 


0.0 


0.1272 


-0.0155 


0.0045 


0.1297 


-0.0091 


0.0095 


0.1277 


05b 


0.4909 


0.0 


0.3727 


0.4941 


0.0187 


0.3701 


0.4909 


0.0095 


0.3723 


06a 


-0.0092 


0.5 


0.1227 


-0.0083 


0.4911 


0.1257 


-0.0081 


0.5046 


0.1246 


06b 


0.4908 


0.5 


0.3772 


0.4842 


0.5042 


0.3730 


0.4919 


0.5046 


0.3754 



structure by about 2.48 eV compared to the cubic phase. 
This result confirms that the monoclinic structure is more 
stable in the low temperature regime. The obtained insu- 
lating gap A g = 0.33 eV is larger than the experimental 
value 0.14 eV, but is close to the values obtained in other 
electronic structure calculations^ One should note that 
the calculated gap is usually larger for the relaxed system 
than for the experimental geometry. 

The charge and spin distribution in the 3d states on 
Fe(B) ions is presented in Table VI. Since we have chosen 
larger Wigner-Seitz radius r-ws = 1-5 A, the obtained 
values are somewhat larger than in other studies . 29 ! 30 
However, the difference between the larger (Bl and BA) 
and smaller (B2 and B3) charges (~ 0.2e) is similar, 
and agrees well with experiment^ The electron density 
of states (DOS) and the orbital density distribution will 
be analyzed in Sec. V, where the EP interactions are 
discussed. 

The parameters of the relaxed crystal structure show 
good agreement with experiment ^ The lattice constants 
are overestimated only by about 1%, and the monoclinic 
angle is well reproduced. As in the real crystal of mag- 
netite, the atoms are displaced predominantly along the 
y direction with respect to the original cubic positions. 
In the monoclinic cell, the y direction corresponds to the 
[110] diagonal one in the Fd3m symmetry. According 
to the experiment^ the largest atomic shift of 0.13 A 
is found for Fe(B3) ions, and it agrees well with the 
theoretical value 0.11 A. Displacements of Fe(B2) and 
Fe(B4) ions are much smaller (0.03-0.04 A), and Fe(-Bl) 
ions do not change their positions at all. As observed 
experimentally, the Fe(^4) ions have only minor displace- 
ments (~ 0.03 A), and they are smaller than those found 
theoretically (0.05 A). The largest shift of oxygen ions 
amounts to 0.11 A, and is slightly larger than the exper- 



imental value found for 04 ions, being 0.09 A. 

This present analysis shows that the largest atomic 
displacements, leading to the monoclinic phase, are con- 
nected with Fe(-B) and O ions. They induce deformations 
of Fes 06 octahedra. It shows that models which assume 
only oxygen vibrations are probably not sufficient to de- 
scribe fully all relevant degrees of freedom involved in the 
Verwey phase transition. In Sees. IV and V we will study 
the phonon spectrum, and in particular we will focus on 
those modes, which strongly influence the electronic and 
crystal structure. 



C. Effect of Hund's exchange 

We close this Section by investigating the effect of 
Hund's exchange on the electronic structure for fixed 
U = 4.0 eV. One finds that the gap value is consider- 
ably enhanced for J = (A g = 0.54 eV) over the value 
found for the realistic exchange interaction J = 0.8 eV 



TABLE VI: Gap in the electronic structure A g (in eV), 
and the charge and magnetization densities for nonequivalent 
Fe(B) atoms, riBi and m,Bi, as obtained in the P2/c phase 
with U — 4.0 eV for two values of Hund's exchange: J = 0.8 
eV (left) and J = (right). 



J 




0.8 






0.0 




A 9 




0.33 






0.54 




atom B 


riBi 




m B i 


UBi 




m B i 


Bl 


6.07 




3.64 


6.09 




3.66 


B2 


5.85 




4.05 


5.82 




4.15 


B3 


5.87 




4.00 


5.85 




4.09 


B4 


6.09 




3.62 


6.11 




3.64 



9 



(A g = 0.54 eV), see Table VI. This demonstrates that the 
gap opens between the minority subbands and the actual 
effective interation between the electrons with the same 
spins is reduced by the exchange term. When Hund's ex- 
change is missing, this interaction is enhanced, resulting 
in an increased gap value. Note that the effective inter- 
action (U — J) between the electrons of the same spin 
decides also about the magnetic and phonon properties 
of the iron metal^ 

The charge density distribution {nsi} over the iron 
ions in B positions is only little modified when J = is 
selected. In fact, the polarization between the amplitude 
of the CO between the {SI, 54} and {B2,BZ} pairs of 
ions is somewhat increased which again confirms that 
the effective interaction between the minority electrons 
has increased. Consequently, the magnetic moments msi 
are also slightly increased for J = over their values 
found for J = 0.8 eV. This behavior allows us to conclude 
that the calculation o 58 ' 79 perfomed with the same value 
of U = 4.0 eV but with larger J = 1.0 eV correspond in 
fact to weaker interaction between the minority electrons. 



IV. LATTICE DYNAMICS 

In this Section, we analyze the lattice dynamics of mag- 
netite in the cubic phase. The low-energy dispersion 
curves were presented already in the previous paper.— 
Here we extend this study by a detailed discussion of 



TABLE VII: Phonon frequencies co n at the zone center com- 
pared with the experimental data of Refs. l38l . l47l . l48l and 
49 (all in meV). I and R denotes infrared and Raman active 
modes, respectively. The difference in phonon frequency due 
to local interactions Slu is given in percents. 



r 


active 


17 = 


U = 4.0 eV 


5oj 


Experiment 






16.84 


17.68 


5.0 


18.5 a 


Ti u 


I 


19.98 


21.46 


7.4 


12 b 


E„ 




21.10 


22.71 


7.0 




T 2g 


R 


24.16 


25.77 


6.7 


23.93 c , 23.93 d 


E s 


R 


32.87 


41.76 


27.0 


39.43 6 , 38.19 c , 37.20 d 


Tx 9 




33.10 


39.27 


18.6 




A 2u 




35.53 


38.08 


7.2 




Tin 


I 


38.31 


40.10 


4.7 


32 6 


Tin 


I 


40.03 


42.85 


7.0 


42. 5 b , 43. 4 C 


Thu 




42.72 


45.31 


6.1 




T 2g 


R 


49.61 


55.22 


11.3 


50.83\ 50.83 d 


E„ 




52.50 


54.30 


3.4 






R 


65.10 


68.89 


5.8 


67.20\ 66.95 c , 66.95 d 


Tin 


I 


66.77 


66.24 


0.8 


68\ 69.43 c 


Ai 9 


R 


73.06 


82.74 


13.2 


83.32 6 , 83.07 c , 82.95 d 


A 2u 




74.21 


81.33 


9.6 





a Reference |3£ 
b Reference 47 
c Reference 4S 
d Reference 4£ 



the phonon modes at the zone center and analyze the ef- 
fects related to electron correlations. In the T point there 
are 42 modes classified according to the IRs of the cubic 
symmetry group, 

r = A lg + 2A 2u + E g + 2E U + T lg + 3T 2g + 5T lu + 2T 2u . 

There are four optic infrared modes with T\ u symmetry, 
five Raman modes with A\ g , T 2g , and E g symmetries, 
respectively, and seven silent modes. In Table VII we 
compare the frequencies of the optic modes with the ex- 
perimental data from the neutron, infrared, and Raman 
measurements. Two sets of theoretical frequencies ob- 
tained from the calculations with: (i) U = and J = 
called below U = 0, and (ii) U = 4.0 eV and J = 0.8 eV 
called below U = 4 eV, are presented. The lowest optic 
mode at 18.5 meV with T 2u symmetry (previously as- 
signed as T 2g ) was measured by the neutron scattering.— 
Four infrared modes with T\ u symmetry have been mea- 
sured by Degiorgi et ai4I Two lowest ones (at 12 and 
32 meV) have very low intensities and large widths in 
experiment, so they can be hardly compared with the 
theoretical values. For the other two modes (at 42.5 and 
67.2 meV) the agreement between the experiment and 
the present calculations is very good. In fact, in a broad 
range the frequencies u) n of these modes are almost inde- 
pendent of U, but finite U is of importance to improve 
the agreement with experiment for the mode observed at 
to = 42.5 meV. 

The theory shows that two frequencies around 40 meV 
are very close to each other, and it may be the reason 
of apparent difficulties in resolving these infrared modes. 
This interpretation could be supported by the fact that 
the peak observed in spectroscopy at 42.5 meV seems 
to be indeed very broad4£ As a success of the present 
theory one finds that all computed Raman modes show 
good agreement with experimental data. In two cases 
(for A\ g and E g modes), calculations with finite U = 4.0 
eV largely improve this agreement. The only discrepancy 
found is the reversed assignement of symmetries (E g and 
T 2g ) of the two modes at ~ 39 and - 50 meV4&22 At 
present we cannot explain this discrepancy and we hope 
that future experiments would verify the present predic- 
tion of the theory. 

Computed phonon frequencies show rather strong de- 
pendence on U . Interestingly, the frequencies increase, in 
spite of larger lattice constant for U > which suggests 
the oposite trend. This indicates a strong effect of lo- 
cal Coulomb interactions on lattice dynamics via respec- 
tive changes in the electronic density distribution and 
in the value of EP coupling. It can be understood us- 
ing a simple picture explained below. In the presence of 
local correlations, the electron localization in Fe(3c?) or- 
bital states is enhanced. It results in weaker screening of 
ionic interactions and, consequently, in larger interatomic 
forces. Since phonon frequencies depend directly on force 
constants, their values increase. This mechanism is sup- 
ported by the observation that the phonon modes stiffen 
in the insulating state in spite of lattice expansion due 
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TABLE VIII: Compatibility relations between the V point 
and A, E and A directions, as well as the X point and A, S 
directions. 



> 



3 19 IL, 



FIG. 2: (Color online) Low-energy phonon frequencies u as 
obtained for the cubic phase of Fe30 4 with: (a) U = J = 
0, and (b) U = 4.0 eV and J = 0.8 eV. The squares show 
the experimental data obtained by neutron scattering.-® Two 
primary OPs are related to A5 and X3 phonons marked by 
circles in (b). The high symmetry points from left to right 
in units of are: L = (3,5,5); T = (0,0,0), X = (0,0,1), 

K = a, ii), r = (1,1,1). 



to the monoclinic distortion i 48 i 49 i 51 One should note also 
that two Raman modes at frequencies ~ 39 meV (E g ) 
and ~ 83 meV (A lg ), which strongly depend on U, show 
the largest anomalies at the transition point i 48 ' 49 In con- 
trast, for the infrared phonons the frequency dependence 
on U is much weaker and no significant anomalies were 
found at TV, apart from some changes in the linewidths 
of these modest 8 - 

The phonon dispersion relations were calculated along 
the main directions of the reciprocal space. In Tab. 
VIII we present the compatibility relations for the IRs in 
the T point and along three directions: A[001], £[110], 
and A[lll]. They show how the degenerate IRs split 
when symmetry is reduced. Knowing these relations, the 
symmetries of all dispersion curves can be properly as- 
signed. Fig. [5] presents the lowest dispersion curves ob- 
tained for U = and U = 4.0 eV, compared with the 
neutron scattering data^ Longitudinal acoustic modes 
agree very well with experimental values, almost inde- 
pendently of the actual value of U (the data for interme- 
diate values of U are not shown). This behavior demon- 
strates that these modes couple very weakly to electronic 



r, x 


A [001] 


E [110] 


A [111] 


Alg 


Ai 


Ei 


Ai 


A 2u 


A 4 


E 3 


Ai 


Eg 
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Xi 
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Ei ©E 3 




x 2 


A 2 © A 3 


E 2 ©E 4 




X, 


A 5 


E 3 ©E 4 




x 4 


A 6 


Ei ©E 2 





density distribution, and are therefore not sensitive to 
changes in electronic structure. A radically different be- 
havior is observed for the transverse acoustic and op- 
tic phonons, which strongly depend on U. The most 
significant changes were obtained here in the A and A 
directions. So far, we could not find a good physical 
explanation for a large discrepancy between theory and 
experiment observed at the L point. 

For the lowest acoustic and optic A 5 modes, the agree- 
ment improves in case of U = 4 eV. It is quite remarkable 
that the optic mode behaves differently in both cases. For 
U = this mode crosses with another A5 mode, so the se- 
quence of phonon branches is interchanged. This behav- 
ior is responsible for the apparent discrepancies between 
the frequencies of phonon modes in theory and experi- 
ment for the U = case. According to the experiment^ 
the lowest phonon at the X point has the X4 symmetry, 
and splits into the Si and £2 phonon branches along the 
[110] direction (see Table VIII). The second lowest mode 
has X3 symmetry and transforms into the £3 and £4 
modes. This order is indeed accurately reproduced by 
the GGA+U calculations with realistic local Coulomb 
interactions (i.e., taking U — 4.0 eV and J = 0.8 eV). 
In contrast, in the U — case, the order of the lowest 
modes is reversed. It shows that Coulomb interactions 
strongly modify lattice dynamics and the phonon energy 
spectrum. We thus emphasize that they have to be in- 
cluded to obtain not only quantitative but even qualita- 
tive agreement with experiment. 



V. ELECTRON-PHONON INTERACTION 

A. Instability of the cubic phase 

Having obtained phonon dispersion curves and their 
symmetries, we now focus on the effects induced by 
phonons, especially on the OPs derived in Sec. II. In 
this Section, we discuss the EP interactions associated 
with the X and A phonons, and present the instabil- 
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ity of the cubic phase triggered by local deformations. 
Since one expects the strongest effects for the low energy 
phonons, we analyze only the lowest energy acoustic and 
optic branches. 

First, for \s.x wave vector we compare the EP coupling 
for the lowest X3 and X4 modes. At the kx point all 
modes are doubly degenerate, so we used polarization 
vectors from only one branch by shifting the wavevector 
to kx + ^ks, removing then the degeneracy. In Fig. [3jwe 
present the total energy E to t of the distorted cubic struc- 
ture as a function of the increasing phonon amplitude Q. 
The energies of relaxed (undistorted) crystals obtained 
in the GGA and GGA+C7 calculations were shifted to 
a common energy origin for a more transparent presen- 
tation (their actual values can be found in Table IV). 
One finds that for U = the energy increases when the 
cubic structure is distorted by either X3 or X4 phonon. 
This behavior is characteristic and represents a typical 
situation when atoms are displaced from their respective 
equilibrium positions. It simply confirms that the cubic 
symmetry is the most stable structure for U — 0, and 
does not undergo any instability towards other possible 
symmetries. 

In the GGA+U calculations, the situation changes 
drastically since the orbital degrees of freedom becomes 
active and influence the electronic structure. As demon- 
strated in Ref. \E§i, the orbital polarization breaks the 
symmetry of wavefunctions within the cubic Fd3m struc- 
ture and leads to lowering of the total energy. This effect 
was demonstrated by initializing the self-consistent pro- 
cedure using the wavefunctions of the C2/C structured 
Applying a similar procedure, we have started calcu- 
lations with the wavefunctions of the low-symmetry 



P2/c phase and optimized the electronic structure self- 
consistently in the Fd3m structure. The reference en- 
ergy obtained by the direct optimization is denoted by 
the point at E to t = in Fig. [3J and the decrease due to 
orbital polarization (about 1.6 eV) is visualised by the 
arrow. The lowering of the ground state energy is as- 
sociated with a partial CO and OO which arises in the 
t2 g states, however, it does not generate the gap opening 
yet. In fact, the calculation suggests that this is rather a 
metastable state, showing only short-range charge-orbital 
correlations without well defined symmetry. 

Total energy E to t decreases further from that obtained 
for the state with CO and OO induced by finite Coulomb 
interactions for increasing phonon distortion Q (see Fig. 
[3j), which leads eventually to the metal-insulator tran- 
sition. The strongest decrease of the total energy was 
found for the X3 mode. On the contrary, by consider- 
ing a single X4 phonon amplitude we have verified that 
this type of distortion does not induce any energy de- 
crease. Instead, the acoustic A5 phonon mode leads to 
the energy lowering when local Coulomb interactions are 
present. Therefore, we included in Fig. [3jalso the energy 
dependence for this acoustic mode. The energy decrease 
is here smaller than in case of the X3 mode. Note how- 
ever that because of double degeneracy the total energy 
may also depend on the phase. Finally, the groundstate 
energy for the P2/c symmetry is shown in Fig. [3] as the 
lowest dashed line. Since the lattice deformation induced 
by either X3 or A5 mode is not sufficient to lower the en- 
ergy to this level, other OPs have to be active as well. 



B. Metal-insulator transition 

The changes in the electron DOS induced by the con- 
sidered phonon modes at the kx point are presented in 
Fig. [U They are compared with the respective DOS 
found for the undistorted structures with FdZm and 
P2/c symmetry, respectively. The results of the calcu- 
lations performed in the GGA and in the GGA+U ap- 
proach are shown on the left and right hand side of Fig. 
21 respectively. The main effect of finite U for the cubic 
structure is to increase the exchange splitting between 
the up and down spin states, cf. Figs. Ufa) and|3Je). 
The up-spin states at Fe(B) ions and down-spin states 
at Fe(.A) ones are shifted to lower energies below Ep. In 
addition, finite U decreases the spectral density of mi- 
nority t2 g states just above Ep. This indicates the en- 
hancement of electron localization due to local electron 
interactions, which leads to gap opening in the distorted 
structure. 

First observation connected with the EP interaction 
is that the overall effects of X3 and X4 distortion are 
much weaker for the uncorrelated (U = 0) case than for 
the correlated (U = 4.0 eV) one. Comparing these two 
distortions for U = [Figs. H^b) and|3Jc)], we see that 
the strongest coupling is associated with the X3 phonon. 
Especially for down-spin states above Ep one finds a de- 
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crease of the spectral density. This instability is con- 
nected with the kx nesting vector at the Fermi surface, 
as discussed some time ago by Yanase and Hamada^i It 
is quite remarkable that this effect alone does not suffice 
to induce a metal-insulator transition, but such a tran- 
sition is triggered when it is amplified by local Hubbard 
interaction U. Indeed, when the X% distortion is made 
in presence of finite U [Fig. [4jf)], the gap opens at Ef, 
and the metal-insulator transition takes place. Note that 
apart from some subtle differences, the distribution of the 
spectral weight near Ep found with the X3 distortion is 
quite similar to that found in the P2/c phase [Fig. ID^h)]. 

Surprisingly, the X4 mode does not produce a similar 
effect, although the changes of the electron DOS were 
here similar to those due to X3 for [7 = 0, and the EP 



coupling is here also much stronger than for [7 = 0. This 
result agrees perfectly with the group theory analysis of 
Sec. II, showing that X3 is the primary OP and X4 
is neither primary nor secondary OP for the VT. We 
remark that the influence of the A5 mode on the DOS has 
been discussed previously^ and it does not induce the 
insulating gap either in spite of its significant coupling 
to the electronic density. Thus, we conclude that the 
X% mode plays the decisive role in the metal-insulator 
transition which accompanies the VT. 
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C. Charge and orbital order 

Realistic treatment of the charge and orbital order re- 
quires the electronic structure calculations with finite lo- 
cal Coulomb interactions using LSDA+[/-like methods. 
The electronic DOS in the cubic phase [Fig. BJe)], as well 
as earlier experimental studies^ suggest that magnetite 
is an itinerant magnet and not a mixed valcnt system. 
The VT is therefore not due to freezing of charge fluc- 
tuations at Fe(-B) ions, but rather due to the symmetry 
change, which leads only to weak CO as a consequence 
of the undergoing phase transition. This scenario is con- 
firmed by the changes in the electron DOS above and 
below the VT, as the gap at the Fermi energy induced 
by the A3 phonon [see Fig. [4jf )] splits the t2 g band into 
the occupied states at the Bl and B2 Fe ions and empty 
states at the B3 and BA ions. Thus, this mode is respon- 
sible for the charge disproportionation which persists also 
in the P2/c phase [Fig. Q|h)]: increased ti g electron den- 
sity for {Bl, B2} ions and reduced density for {B3, B4} 
ions. Note that four Fe(-B) sites split into two subclasses 
as a result of the lattice distortion in the A3 mode [Fig. 
Eta)]. The amplitude of the charge (spin) disproportion- 
ation is comparable with that found for the P2/c sym- 
metry (see Table VI). Interestingly, this CO pattern is 
similar to that proposed originally by Verwey and fulfills 
the Anderson criterion. These findings agree also with 
other LSDA+C/ and GGA+C/ calculations-^^!™ 

In addition, we have found that the A3 mode stabilizes 
the 00 in the ti g states. The resulting orbital polariza- 
tion with atomic displacements is presented in Fig. EJa). 
In the considered mode the atoms vibrate along the [110] 
direction, and the dominating components are those con- 
nected with distortions of Fe(B) ions and oxygen ions in 
the same planes. Oxygen displacements are 65% of Fc(B) 
displacements, while Fe(A) atoms shift only by about 
23% of them (not shown). The considered mode is lim- 
ited to single Fe-0 planes: Fe(B) ions move along with O 
ions in one plane, while atoms in neighboring planes do 
not move at all (they participate in the other A3 branch 
in the perpendicular direction). This deformation modi- 
fies distances between Fe(i?) atoms and oxygens situated 
above and below them, inducing changes in the electron 
density distribution within t2 g states and a coexisting 
charge and orbital order in the Bl and B2 chains. As 
indicated in Fig. [5] electrons occupy orthogonal d yz and 
d xz orbitals, thus forming the state with alternating 00. 

This coexisting charge-orbital order and the insulat- 
ing gap induced by the A3 mode has to be compared 
with that observed in the monoclinic phase of P2 / c sym- 
metry [Fig. EJb)]. The electronic DOS for this sym- 
metry is presented in Fig. H^h). The magnitude of the 
gap is similar to the gap generated by the A3 phonon 
mode shown in Fig. 0Jf), however, the occupied states 
below Ep belong now to the {B1,BA} ions and empty 
states to the {B2,B3} ones. Actually, Fe(B) ions are 
split into two groups {Bl, BA} and {B2, B3}, as average 
Fe-0 distances for the bonds to the {B2,B3} sites are 
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significantly smaller than those for the {B1,B4} ones. 
Therefore, the CO has to change from that promoted by 
the A3 mode alone. Surprisingly, the weak CO found in 
the LT monoclinic phase does not satisfy the Anderson 
criterion, but plays no role in the actual mechanism of 
the transition. 

The corresponding orbital pattern also has changed 
[see Fig. Hp 3 )]. The occupied states on Bl and BA 
irons form the alternating 00 with the orthogonal or- 
bitals on the nearest neighbor iron sites. We emphasize 
that in the P2/c frame, the occupied ti g states are given 
by -^(d xz ±d y z) and d x 2_ y 2 combinations. Therefore, on 

comparing Figs. [5ja) andEtb), we conclude that the A 3 
mode demonstrates a generic tendency towards the 00 
within the t2 g minority states, but only partly explains 
the observed 00 in the LT phase with P2/c symmetry. 
This is however quite natural and expected, since the A3 
mode is only one component of the OPs that condense at 
the VT, which all finally lead to the change of symmetry 
to the LT monoclinic phase. 

In the next step, one should consider the superposi- 
tion of the A3 with the A5 mode. The latter one is 
crucial to explain the doubling of the unit cell and the 
observed charge-orbital pattern. As we mentioned before, 
this mode is double degenerate and the atomic displace- 
ments depend on the phase. So, a detailed analysis of the 
EP coupling, similar to that presented above for the A3 
modes, is rather difficult. Moreover, because of the sys- 
tem complexity, it is likely that the atomic displacements 
belonging to the secondary OPs are also involved in the 
final charge-orbital order. In spite of these difficulties, 
the studies presented in this Section demonstrate how 
the phonons interplay with charge and orbital degrees of 
freedom in the VT. 
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VI. MECHANISM OF THE VERWEY 
TRANSITION 

From the results presented in the previous Section and 
from the earlier experimental studies, it is clear that the 
mechanism of the VT includes two essential ingredients: 
(i) strong intraatomic Coulomb interactions {U, J} at Fe 
ions, and (ii) phonon-driven lattice instability. On the 
one hand, the electron interactions are responsible for 
the orbital polarization, which breaks the symmetry of 
the wavefunction and enhances the electron localization 
in the <2g states. On the other hand, the phonon OPs 
induce the crystal distortion, which generates the struc- 
tural transformation from the cubic Fd3m phase to the 
LT phase with monoclinic P2/c symmetry. The coupling 
between these two subsystems is crucial since none of 
them alone would be able to explain the metal-insulator 
transition which occurs in magnetite at Ty. As we have 
shown, without the Hubbard interaction U, the EP inter- 
action would be too weak to induce structural changes. 
In the presence of electron correlations, the coupling be- 
tween the electrons in t2 g states and phonons is largely 
enhanced since it leads to stabilization of the 00 and to 
lowering of the total energy. The corresponding CO is 
very subtle and seems to be rather a consequence of the 
joint effect of the 00 and enhanced EP interactions than 
the driving force of the VT. 

The splitting of the t2 g states due to the EP interac- 
tion is usually referred to as the Jahn- Teller effect. In 
principle, the mechanism observed in magnetite is simi- 
lar to this effect, however, it differs in many respects from 
the classical examples. In a typical situation, one consid- 
ers the displacements of oxygen ions around the transi- 
tion metal cations, which break the symmetry and split 
the degenerate 3d states. As we discussed in the previ- 
ous Sections, the displacement pattern in the monoclinic 
structure is more complicated and involves also Fe(B) 
ions. In particular, in the A3 mode the iron atoms (with 
largest amplitudes) vibrate in the same plane as oxy- 
gens. Therefore, the orbital polarization results in this 
case mainly from the interplanar interactions between 
iron and oxygen ions. 

A dynamical coupling between the charge-orbital fluc- 
tuations and phonons induces the critical diffuse scatter- 
ing above Ty . Temperature dependence of this scattering 
gives us important information about precursors, which 
show up far from the critical point. The signs of the tran- 
sition appear already at about 200 K. Apart from neutron 
scattering coming from kr and kx points ; 43 ' 45 there are 
maxima at other positions, also with an incommensurate 
wavelengths. For instance, very pronounced critical be- 
havior was observed at k = (8, 0, 0.75)42 When temper- 
ature is lowered, the coupling between different modes, 
described by the Landau free energy (|T2"]). leads to stabi- 
lization of particular phonons with commensurate wave- 
lengths. Very close to the transition point, about 5 K 
above Ty, the A 5 OP becomes active, producing the sig- 
nal in neutron measurement* 3 ^ The question arises why 



is this mode observed only a few kelvins above Ty , while 
signals from the T and X points are seen at much higher 
temperatures? This may be connected with the relative 
strength of different OPs. In Fig. [21 we see that the 
coupling to A5 mode is weaker than to A3 mode, and it 
gives a different signal in neutron scattering. 

Critical scattering is one of indications of a short-range 
(charge and orbital) order above the VT, first discussed 
by Andersoni 3 - Another evidence is connected with the 
change in entropy, which is smaller than in the order- 
disorder phase transitions. 86 Also the temperature de- 
pendence of resistivity at T > Ty, which is not typical 
for metals, suggests the polaronic nature of charge carri- 
ers. First photoemision measurements indicated the gap 
closing above Ty^ but more detailed studie a 33 ' 34 ' 87 es- 
tablished that the gap only decreases, without any sharp 
change at Ty. The optical conductivity studies suggest 
also the gap opening below Ty, but above the transi- 
tion the conductivity spectrum does not exhibit a metal- 
lic Drude-type behavior. 56 Instead, a hopping type con- 
ductivity is observed, with highly diffuse character of 
charge dynamics. This picture is supported by the struc- 
tural EXAFS studies, which found that the local crys- 
tal geometry does not change at Ty, 50 It means that 
the quasistatic lattice distortions are present already in 
the high-symmetry cubic phase above Ty. All these re- 
sults are consistent with the present theory, which shows 
that some phonons strongly couple to electronic states 
and may induce local crystal deformations and polaronic 
short range order above Ty. 

As demonstrated in this work and the previous study^ 
the theoretical approach based on the GGA+[/ (or 
LDA+Z7) provides a very good description of the VT. 
However, there are some aspects, which involve dynami- 
cal processes and many-body interactions. For example, 
the detailed analysis of the photoemision spectra showed 
that the DFT is not sufficient to explain changes in- 
duced by the VT and better agreement with experiment 
can be achieved only by using the DMFT approach. 58 
An interesting effect was observed in the spin excitation 
spectrum, which shows a large splitting in the acoustic 
magnon branch at kA below the VT. 88 This effect cannot 
be explained merely by the dependence of exchange inter- 
actions on the crystal structure or charge ordering, but 
is rather a consequence of the magnon-phonon coupling. 
Additionally, a recent neutron scattering study revealed 
an anomalous broadening and energy shift of the A5 spin 
wave above Ty^ A strong spin-phonon interaction, sug- 
gested by the present work, provides a sound starting 
point for studying such effects. 

Finally, the behavior of magnetite at high pressure is 
rather intriguing. According to diffraction studies^ the 
critical temperature Ty decreases with pressure, so the 
VT can be induced at rather low temperature by apply- 
ing pressure in the LT monoclinic phase. Recently two 
contradicting views on the VT based on pressure experi- 
ments have been suggested. On the one hand, a possible 
transition from the inverse to normal spinel structure has 
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been suggested as consistent with the interpretation of 
the Mossbauer and x-ray scattering measurements ! 91 ! 92 
This scenario precludes a CO on the Fe(B) sites below 
TV, but predicts a large ~ 50% increase in the bulk mag- 
netization with increasing pressure or by lowering tem- 
perature. On the other hand, no change in magnetic 
moments at both A and B sites has been detected by the 
neutron diffraction studies in a broad range of pressure 
up to p = 5.3 GPa^ This latter result rules out a posibil- 
ity of the inverse spinel to normal spinel transition under 
decreasing temperature (increasing pressure), at least in 
the regime of pressure lower than 5.3 GPa. It also agrees 
with earlier band structure calculations and with the 
present mechanism of structural transition in magnetite. 
As the electron occupations at A and B sites are modi- 
fied rather weakly, no significant change in the values of 
magnetic local moments is expected. We emphasize that 
this latter view does not imply large CO below the VT. 
In fact, the structural transition suggested in the present 
work not only explains the dramatic change in the con- 
ductivity observed at Ty, but also predicts rather weak 
CO, in agreement with experiment j 18 i 19 ' 20 i 21 

VII. SUMMARY AND CONCLUSIONS 

In this work, we have presented a detailed group the- 
ory analysis of the Verwey transition. We have identified 
three primary OPs: X3, A5, and Tig, which describe 
the symmetry reduction in the crystal structure trans- 
formation from the high-temperature cubic Fd3m phase 
to the LT monoclinic P2/c phase. By performing the 
numerical ab initio computations, we have demonstrated 
that a prominent role is played by the A3 mode which: 
(i) couples strongly to the electronic states, (ii) lowers 
the total energy, and (iii) is responsible for the metal- 
insulator transition. The latter transition occurs only 
when local Coulomb interaction U is explicitly included 
in tig iron states and may trigger weak charge order at B 
sites. In fact, this pnenomenon appears counterintuitive 
and occurs only as a result of the accompanying local 
lattice distortions which may be seen as a Jahn- Teller 
lattice instability. These results support the recent point 
of vie w 60 ' 94 that neither pure electrostatics is the main 
factor responsible for the charge order observed below the 
Verwey transition, nor the charge order is the mechanism 
driving the transition by itself. At the same time, the 
local Coulomb interaction generates alternating orbital 
order, which leads to strong reduction of charge mobility 



and amplifies electron-lattice effects. 

The present study reconciles several previous points of 
view on the Verwey transition in magnetite and suggests 
that the physical effects which occur simultaneously be- 
low this transition in the monoclinic P2/c phase can be 
classified into the ones which are the primary cause of 
the symmetry change and the ones which occur only as 
its consequence. In this way it contributes to the recent 
debate concerning the origin of the transition and clari- 
fies the role played in it by the charge order. While weak 
charge order has been found at Fe(B) sites, it is not sur- 
prising that it does not obey the Verwey model. In fact, 
it is only one of the manifestations of strong local electron 
interactions in partly occupied ti g states rather than the 
primary cause of the observed symmetry change. 

We have also compared phonon energy spectrum of 
magnetite with the experimental data obtained by Ra- 
man and infrared spectroscopy, as well as by neutron 
scattering. We have found that phonon frequencies 
strongly depend on local electron interactions, and the 
GGA+U calculations performed with realistic parame- 
ters for Fe ions (U — 4.0 eV and J = 0.8 eV) give very 
satisfactory qualitative and quantitative agreement with 
the experimental data. In contrast, when the electron 
interaction effects are neglected, the phonon spectra are 
even qualitatively different from the observed ones. 

Summarizing, we have shown that the Verwey transi- 
tion is promoted by a set of order parameters with mixed 
electron and phonon character. The insulating mono- 
clinic P2 1 c phase occurs below the transition as a result 
of the instability driven by the electron-phonon coupling 
in presence of strong electron correlations. We argue that 
the electron-lattice coupling plays an important role also 
in other transition metal oxides with strongly correlated 
electrons. Further studies of the lattice relaxation effects 
in these systems may lead to discoveries of new electronic 
phenomena that could be understood only by simultane- 
ous treatment of electronic and lattice degrees of freedom. 
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